
# Description -------------------------------------------------------------


# Replication Code for the Paper 

# "Rebels, Revenue, and Redistribution:
# The Political Geography of Post-Conflict
# Power-Sharing in Africa"

# Authors: Felix Haass, Martin Ottmann



# Libraries ---------------------------------------------------------------

library(lfe)
library(sf)
library(tidyverse)
library(stargazer)
library(countrycode)
library(haven)
library(caret)
library(doParallel)

# Custom Functions ---------------------------------------------------------

# predict values from fixed effects linear models
source("./code/custom_functions/predict_felm.R")

# Load Data ---------------------------------------------------------------

# this loads the data frame "rrr_data", the main data frame
load("./data/rrr_data.rda")

# load PSED data 
load("./data/psed_month.rda")

#  ** TABLES ** -----------------------------------------------------------


# Table 2 -----------------------------------------------------------------


base_fe_notime <- felm(log(nlights_1years ) ~
                         cabinetINC 
                       
                       
                       | gid  | 0 | gid , 
                       data = rrr_data , keepCX = T)


base_fe_time <- felm(log(nlights_1years ) ~
                       cabinetINC 
                     
                     
                     | gid_fct + country_year | 0 | gid , 
                     data = rrr_data, keepCX = T)


base_fe_timecontrols <- felm(log(nlights_1years) ~
                               cabinetINC  + 
                               
                               log(pop_ip + 1) +
                               log(brd_rollsum + 1)  +
                               log(nonstate_rollsum +1 )
                             
                             | gid_fct + country_year  | 0 | gid  , 
                             data = rrr_data, keepX = T)


base_fe_time_controlslagfe <- felm(log(nlights_1years ) ~
                                     cabinetINC +
                                     log(pop_ip + 1) +
                                     log(brd_rollsum + 1)  +
                                     log(nonstate_rollsum +1 ) +
                                     log(nlights_tmin1 + 0.01)  
                                   
                                   
                                   | gid + country_year | 0 | gid, 
                                   data = rrr_data , keepCX = T)

base_fe_timecontrols_inc_control <- felm(log(nlights_1years) ~
                                           cabinetINC  + 
                                           
                                           log(pop_ip + 1) +
                                           log(brd_rollsum + 1)  +
                                           log(nonstate_rollsum +1 )
                                         
                                         | gid_fct + country_year | 0 | gid, 
                                         data = rrr_data %>% 
                                           group_by(gid) %>% 
                                           filter(sum(eprINC) == 6 | cabinetINC_treatgroup == 1),  
                                         keepX = T)

base_fe_timecontrols_exc_control <- felm(log(nlights_1years) ~
                                           cabinetINC  + 
                                           
                                           log(pop_ip + 1) +
                                           log(brd_rollsum + 1)  +
                                           log(nonstate_rollsum +1 )
                                         
                                         | gid_fct + country_year | 0 | gid, 
                                         data = rrr_data %>%
                                           group_by(gid) %>% 
                                           filter(sum(eprINC) == 0 | cabinetINC_treatgroup == 1),  
                                         keepX = T)

# * Output  -----------------------------------------------------------

stargazer(base_fe_notime, 
          base_fe_time,
          base_fe_timecontrols, 
          base_fe_time_controlslagfe,
          base_fe_timecontrols_exc_control,
          base_fe_timecontrols_inc_control, 
          type = "html",
          out = "./tables/table2.html", 
          
          covariate.labels = c("Representation in Power-Sharing", 
                               "Population (log)", 
                               "Past Battle Fatalities (log)", 
                               "Past Non-State Fatalities (log)", 
                               "Night Lights$_{t-1}$ (log)"),
          keep.stat = c("n", "adj.rsq"),
          add.lines = list(c("Grid-cell FE", 
                             rep("\\multicolumn{1}{c}{Yes}", 6)), 
                           c("Country-Year FE", 
                             "\\multicolumn{1}{c}{-}",
                             rep("\\multicolumn{1}{c}{Yes}", 5)), 
                           c("Control Group", 
                             "\\multicolumn{1}{c}{All}", 
                             "\\multicolumn{1}{c}{All}", 
                             "\\multicolumn{1}{c}{All}",
                             "\\multicolumn{1}{c}{All}",
                             "\\multicolumn{1}{c}{Excluded}",
                             "\\multicolumn{1}{c}{Included}")),
          float = F,
          align = T, 
          notes.align = "l",
          dep.var.labels = "Night Lights$_{t+1}$ (log)",
          dep.var.caption = NULL, 
          notes        = "", 
          notes.label = "", 
          notes.append = FALSE, 
          out.header = F)


# ** FIGURES ** -----------------------------------------------------------


# Figure 1 ----------------------------------------------------------------

# some necessary data transformation from raw PSED data
psed_cy <- psed_month %>% 
  group_by(location, pcm1) %>% 
  summarise(year_had_cabINC = ifelse(1 %in% cabinetINC, 1, 0), 
            year_had_tpsINC = ifelse(1 %in% territorialINC, 1, 0),
            year_had_milINC = ifelse(1 %in% milintegrationEVE, 1, 0), 
            year_had_epsINC = ifelse(1 %in% epscommissionINC | 1 %in% epscompanyINC, 1, 0)) %>% 
  mutate(region = countrycode(location, "country.name", "continent"))

ps_types_in_africa_per_year <- psed_cy %>% 
  group_by(region, pcm1) %>% 
  summarise(instances_of_cabPS = sum(year_had_cabINC), 
            instances_of_terrPS = sum(year_had_tpsINC), 
            instances_of_milINC = sum(year_had_milINC), 
            instances_of_eps = sum(year_had_epsINC)) %>% 
  gather(key = type_of_ps, value = power_sharing, -pcm1, -region) %>% 
  rename(pcm1 = pcm1) 

# generate plot
plot_output <- ps_types_in_africa_per_year %>% 
  filter(region != "Oceania") %>% 
  ggplot(. , 
         aes(pcm1, 
             y = power_sharing, 
             fill = region,
             group = region)) + 
  geom_area(color = "black", size = 0.3) +
  geom_line(size = 0.3) +
  scale_fill_grey() + 
  scale_y_continuous(breaks = seq(0, 10, 2)) +
  facet_grid(region ~ type_of_ps, 
             labeller = labeller(.cols = c(instances_of_cabPS = "Executive Power-Sharing", 
                                           instances_of_terrPS = "Territorial Power-Sharing", 
                                           instances_of_milINC = "Military Power-Sharing", 
                                           instances_of_eps = "Economic Power-Sharing"))) +
  theme_bw() +
  theme(legend.position = "none") +
  theme(strip.background = element_rect(fill = "lightgrey")) +
  theme(text = element_text(size=8), 
        panel.grid = element_blank()) +
  labs(y = "Instances of Power-Sharing Practices", 
       x = "Post-Conflict Month")


# * Save Figure -----------------------------------------------------------

ggsave(plot_output, file = "./plots/figure1.pdf", 
       width = 8, height = 5, scale = 0.8)



# Figure 2 (a) ------------------------------------------------------------


# * Load Africa Data ---------------------------------------------------

load("./data/map_data/africa_map.rda")
load("./data/map_data/priogrid_sample_countries.rda")
load("./data/map_data/rrr_countries.rda")


# * Generate Plot ---------------------------------------------------------

africa_map_rrr <- ggplot() +
  geom_sf(data = africa_map %>% st_as_sf(), fill = NA, color = "grey", size = 0.3) +
  geom_sf(data = priogrid_samplecountries %>% st_as_sf(), alpha = 0.9, 
          aes(fill = as.factor(cabinetINC_treatgroup)), 
          size = 0.08) +
  geom_sf(data = rrr_countries %>% st_as_sf(), fill = NA, color = "black", size = 0.4) +
  coord_sf(datum = NA) +
  scale_fill_manual("Rebels' ethnic \nconstituency area", labels = c("No", "Yes"), 
                    values = c("#fee0d2", "#de2d26")) + 
  guides(fill = guide_legend(title.position="top")) +
  theme_minimal() +
  theme(legend.position = "bottom", 
        legend.justification = "center") 


# * Save Plot -------------------------------------------------------------


ggsave(africa_map_rrr, filename = "./plots/figure2a.pdf", 
       height = 7.5, width = 8)


# Figure 2 (b) ------------------------------------------------------------

# * Load NL data

load("./data/map_data/priogrid_nlights.rda")


# * Generate Plot ---------------------------------------------------------

africa_map_nl <- ggplot() +
  geom_sf(data = africa_map %>% st_as_sf(), fill = NA, color = "grey", size = 0.3) +
  
  geom_sf(data = priogrid_nlights %>% st_as_sf(), 
          aes(fill = nl_change > 0), 
          size = 0.08) +
  geom_sf(data = rrr_countries %>% st_as_sf(), fill = NA, color = "black", size = 0.4) +
  
  scale_fill_manual("Night lights two years \nafter start of power-sharing", 
                    values = c("Black", "White"), 
                    labels = c("No increase", "Increase")) +
  guides(fill = guide_legend(title.position="top")) +
  coord_sf(datum = NA) +
  theme_minimal() +
  theme(legend.position = "bottom", 
        legend.justification = "center") 


# * Save Plot -------------------------------------------------------------

ggsave(africa_map_nl, filename = "./plots/figure2b.pdf", 
       height = 7.5, width = 8)



# Figure 3: Data Illustration ---------------------------------------------

# Fully replicating figure 3 requires loading of about 10 raw night light files 
# totalling about 6.7GB of satellite data. Moreover, a custom Google API key is
# necessary to download the background map data. To keep the replication file 
# size manageable, we do not include the data + code in this file. 

# We are happy to provide the code + raw NL data upon request: haass.felix@gmail.com

# Figure 4: Cross Validation---------------------------------------------------


# * Prepare data for prediction ---------------------------------------------

rrr_ml_data <- rrr_data %>% 
  select(gid, gid_fct, nlights_1years, cabinetINC, cabinetINC_leader,
         pop_ip, brd_rollsum, nonstate_rollsum, nlights_tmin1, 
         petroleum_dummy_cf, nlights_calib_mean,
         urban_cf, agri_gc, capdist, 
         country_year) %>% 
  filter(!is.na(urban_cf)) %>% 
  # log transform variables
  mutate(nlights_1years = log(nlights_1years), 
         pop_ip = log(pop_ip +1), 
         brd_rollsum = log(brd_rollsum + 1), 
         capdist = log(capdist + 1), 
         urban_cf = log(urban_cf + 0.1), 
         nonstate_rollsum = log(nonstate_rollsum + 1), 
         nlights_tmin1 = log(nlights_calib_mean + 0.01))

# we do a group K Fold sampling on grid cells, in order to split by grid cell 
# so the same grid is  not used both in test and training data

# this generates the holdouts
gid_kfold <- groupKFold(rrr_ml_data$gid, k = 15)

tc <- trainControl(
  method = "cv", # cross validation
  number = 10,
  verboseIter = T, 
  allowParallel = T, 
  index = gid_kfold
)

# use multi core 
cl <- makePSOCKcluster(4)

registerDoParallel(cl)

# specify formula
formula_full <- nlights_1years ~ cabinetINC + 
  pop_ip +
  nlights_tmin1 + 
  agri_gc + 
  urban_cf + 
  capdist + 
  brd_rollsum + 
  nonstate_rollsum +
  petroleum_dummy_cf

method_string <- "lm"

# estimate model
model_full <- train(formula_full, 
                    data = rrr_ml_data, 
                    trControl = tc, 
                    method = method_string)

# stop multi core
stopCluster(cl)


model_full_importance <- varImp(model_full, scale = F) 

# * Extract Variable importance -------------------------------------------

model_full_importance <- data.frame(model_full_importance$importance) %>% 
  mutate(variables = row.names(.)) %>% 
  mutate(variables = forcats::fct_reorder(variables, Overall), 
         start = 0)


# * Generate Plot ---------------------------------------------------------

variable_importance <- ggplot(model_full_importance, 
                              aes(x = variables, y = Overall)) +
  geom_point(size = 1.3) +
  geom_segment(aes(x=variables, 
                   xend=variables, 
                   y=0, 
                   yend=Overall)) +
  scale_x_discrete(labels = c("Battle Fatalities", 
                              "Non-State Fatalities", 
                              "Distance to Capital", 
                              "Petroleum Deposit", 
                              "Land Use (Agriculture)", 
                              "Urbanization", 
                              "Executive Power-Sharing", 
                              "Population", 
                              "Night Lights (t-1)")) +
  labs(x = "", y = "Variable Importance") +
  theme_bw()+
  theme(panel.grid = element_blank()) +
  coord_flip()



# * Save Plot -------------------------------------------------------------


ggsave(variable_importance, filename = "./plots/figure4.pdf", 
       height = 8, width = 6, scale = 0.8)



# Figure 5 ----------------------------------------------------------------


# * Any Cabinet Level------------------------------------------------------

base <- felm(log(nlights_calib_mean + 0.01) ~
               cabinetINC + 
               
               log(pop_ip + 1) +
               log(brd_rollsum + 1)  +
               log(nonstate_rollsum +1 )
             
             
             | gid_fct + country_year  | 0 | gid, 
             data = rrr_data,keepCX = T)

base_1years <- felm(log(nlights_1years) ~
                      cabinetINC + 
                      
                      log(pop_ip + 1) +
                      log(brd_rollsum + 1)  +
                      log(nonstate_rollsum +1 )
                    
                    
                    | gid_fct + country_year  | 0 | gid, 
                    data = rrr_data,keepCX = T)


base_2years <- felm(log(nlights_2years) ~
                      cabinetINC + 
                      
                      log(pop_ip + 1) +
                      log(brd_rollsum + 1)  +
                      log(nonstate_rollsum +1 )
                    
                    
                    | gid_fct + country_year  | 0 | gid, 
                    data = rrr_data,keepCX = T)

base_3years <- felm(log(nlights_3years) ~
                      cabinetINC + 
                      
                      log(pop_ip + 1) +
                      log(brd_rollsum + 1)  +
                      log(nonstate_rollsum +1 )
                    
                    
                    | gid_fct + country_year  | 0 | gid, 
                    data = rrr_data,keepCX = T)

base <- broom::tidy(base) %>% 
  mutate(year = "+/- 0", type = "cabinetINC")

base_1years <- broom::tidy(base_1years) %>% 
  mutate(year = "+ 1", type = "cabinetINC")

base_2years <- broom::tidy(base_2years) %>% 
  mutate(year = "+ 2", type = "cabinetINC")

base_3years <- broom::tidy(base_3years) %>% 
  mutate(year = "+ 3", type = "cabinetINC")

all_time_models_cabINC <- bind_rows(base, base_1years, base_2years, base_3years) %>% 
  filter(term == "cabinetINC")


# * Econ Portfolio---------------------------------------------------



base <- felm(log(nlights_calib_mean + 0.01) ~
               econcabINC + 
               
               log(pop_ip + 1) +
               log(brd_rollsum + 1)  +
               log(nonstate_rollsum +1 )
             
             
             | gid_fct + country_year  | 0 | gid, 
             data = rrr_data,keepCX = T)

base_1years <- felm(log(nlights_1years) ~
                      econcabINC + 
                      
                      log(pop_ip + 1) +
                      log(brd_rollsum + 1)  +
                      log(nonstate_rollsum +1 )
                    
                    
                    | gid_fct + country_year  | 0 | gid, 
                    data = rrr_data,keepCX = T)


base_2years <- felm(log(nlights_2years) ~
                      econcabINC + 
                      
                      log(pop_ip + 1) +
                      log(brd_rollsum + 1)  +
                      log(nonstate_rollsum +1 )
                    
                    
                    | gid_fct + country_year  | 0 | gid, 
                    data = rrr_data,keepCX = T)

base_3years <- felm(log(nlights_3years) ~
                      econcabINC + 
                      
                      
                      
                      log(pop_ip + 1) +
                      log(brd_rollsum + 1)  +
                      log(nonstate_rollsum +1 )
                    
                    
                    | gid_fct + country_year  | 0 | gid, 
                    data = rrr_data,keepCX = T)

base <- broom::tidy(base) %>% 
  mutate(year = "+/- 0",type = "econcabINC")

base_1years <- broom::tidy(base_1years) %>% 
  mutate(year = "+ 1", type = "econcabINC")

base_2years <- broom::tidy(base_2years) %>% 
  mutate(year = "+ 2", type = "econcabINC")

base_3years <- broom::tidy(base_3years) %>% 
  mutate(year = "+ 3", type = "econcabINC")

all_time_models_econINC <- bind_rows(base, base_1years, base_2years, base_3years) %>% 
  filter(term == "econcabINC")


# * Senior Portfolio ---------------------------------------------------------------



base <- felm(log(nlights_calib_mean  + 0.01) ~
               seniorINC + 
               
               
               
               log(pop_ip + 1) +
               log(brd_rollsum + 1)  +
               log(nonstate_rollsum +1 )
             
             
             | gid_fct + country_year  | 0 | gid, 
             data = rrr_data,keepCX = T)

base_1years <- felm(log(nlights_1years) ~
                      seniorINC + 
                      
                      
                      
                      log(pop_ip + 1) +
                      log(brd_rollsum + 1)  +
                      log(nonstate_rollsum +1 )
                    
                    
                    | gid_fct + country_year  | 0 | gid, 
                    data = rrr_data,keepCX = T)


base_2years <- felm(log(nlights_2years) ~
                      seniorINC + 
                      
                      
                      
                      log(pop_ip + 1) +
                      log(brd_rollsum + 1)  +
                      log(nonstate_rollsum +1 )
                    
                    
                    | gid_fct + country_year  | 0 | gid, 
                    data = rrr_data,keepCX = T)

base_3years <- felm(log(nlights_3years) ~
                      seniorINC + 
                      
                      
                      
                      log(pop_ip + 1) +
                      log(brd_rollsum + 1)  +
                      log(nonstate_rollsum +1 )
                    
                    
                    | gid_fct + country_year  | 0 | gid, 
                    data = rrr_data,keepCX = T)

base <- broom::tidy(base) %>% 
  mutate(year = "+/- 0",type = "seniorINC")

base_1years <- broom::tidy(base_1years) %>% 
  mutate(year = "+ 1", type = "seniorINC")

base_2years <- broom::tidy(base_2years) %>% 
  mutate(year = "+ 2", type = "seniorINC")

base_3years <- broom::tidy(base_3years) %>% 
  mutate(year = "+ 3", type = "seniorINC")

all_time_models_seniorINC <- bind_rows(base, base_1years, base_2years, base_3years) %>% 
  filter(term == "seniorINC")



# * Leader Constituency-----------------------------------------------



base <- felm(log(nlights_calib_mean + 0.01 ) ~
               cabinetINC_leader + 
               
               
               
               log(pop_ip + 1) +
               log(brd_rollsum + 1)  +
               log(nonstate_rollsum +1 )
             
             
             | gid_fct + country_year  | 0 | gid, 
             data = rrr_data,keepCX = T)

base_1years <- felm(log(nlights_1years) ~
                      cabinetINC_leader + 
                      
                      
                      
                      log(pop_ip + 1) +
                      log(brd_rollsum + 1)  +
                      log(nonstate_rollsum +1 )
                    
                    
                    | gid_fct + country_year  | 0 | gid, 
                    data = rrr_data,keepCX = T)


base_2years <- felm(log(nlights_2years) ~
                      cabinetINC_leader + 
                      
                      
                      
                      log(pop_ip + 1) +
                      log(brd_rollsum + 1)  +
                      log(nonstate_rollsum +1 )
                    
                    
                    | gid_fct + country_year  | 0 | gid, 
                    data = rrr_data,keepCX = T)

base_3years <- felm(log(nlights_3years) ~
                      cabinetINC_leader + 
                      
                      
                      
                      log(pop_ip + 1) +
                      log(brd_rollsum + 1)  +
                      log(nonstate_rollsum +1 )
                    
                    
                    | gid_fct + country_year  | 0 | gid, 
                    data = rrr_data,keepCX = T)

base <- broom::tidy(base) %>% 
  mutate(year = "+/- 0", type = "cabinetINC_leader")

base_1years <- broom::tidy(base_1years) %>% 
  mutate(year = "+ 1", type = "cabinetINC_leader")

base_2years <- broom::tidy(base_2years) %>% 
  mutate(year = "+ 2", type = "cabinetINC_leader")

base_3years <- broom::tidy(base_3years) %>% 
  mutate(year = "+ 3", type = "cabinetINC_leader")

all_time_models_cabinetINC_leader <- bind_rows(base, base_1years, base_2years, base_3years) %>% 
  filter(term == "cabinetINC_leader")




# Extract Model Data for Plotting ---------------------------------------------


all_models <- bind_rows(all_time_models_cabINC, 
                        all_time_models_cabinetINC_leader, 
                        all_time_models_econINC, 
                        all_time_models_seniorINC) %>% 
  mutate(year = forcats::fct_relevel(year, c("+/- 0", "+ 1", "+ 2", "+ 3"))) %>% 
  mutate(type = forcats::fct_relevel(type, c("cabinetINC", "seniorINC", "econcabINC", "cabinetINC_leader")))



# * Generate Plot ---------------------------------------------------------

rrr_viz_results <- ggplot(all_models, 
                          aes(x = year, group = type)) + 
  geom_point(aes(y = estimate), 
             position = position_dodge(0.5)) +
  geom_errorbar(aes(ymin = estimate - 1.96 * std.error, 
                    ymax = estimate + 1.96 * std.error), 
                width = 0,
                position = position_dodge(0.5)) +
  facet_grid(~type, labeller = labeller(type = c(cabinetINC = "Any Cabinet Level", 
                                                 econcabINC = "Economic Portfolio", 
                                                 
                                                 cabinetINC_leader = "Leader Constituency", 
                                                 seniorINC = "Senior Portfolio"))) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  labs(x = "Years after start of power-sharing government", y = "Coefficient Estimate") +
  theme_bw() +
  theme(panel.grid = element_blank(),
        legend.position = "none")


# * Save Plot ---------------------------------------------------------------

ggsave(rrr_viz_results, filename = "./plots/figure5.pdf", 
       width = 9, height = 4, scale = .9)


# Figure 6 ----------------------------------------------------------------


# estimate interaction model
model_group_interaction <- felm(log(nlights_1years) ~
                                  cabinetINC + 
                                  cabinetINC:number_of_groups_na_rm 
                                
                                
                                | gid + country_year  | 0 | gid, 
                                data = rrr_data, keepX = T)

# predict values
predictions_inter <- predict.felm(model_group_interaction, 
                                  newdata = data.frame(cabinetINC = c(rep(1, 5)), 
                                                       number_of_groups_na_rm = c(1:5)), 
                                  se.fit = T)

# data for plotting
predictions_inter <- data.frame(estimate = predictions_inter$fit, 
                                ub = predictions_inter$fit + 1.96 * predictions_inter$se.fit, 
                                lb = predictions_inter$fit - 1.96 * predictions_inter$se.fit, 
                                number_of_groups = 1:5)



# * Generate Plot ---------------------------------------------------------

num_groups_plot <- ggplot(predictions_inter, aes(x = number_of_groups, y = fit)) + 
  geom_point() + 
  geom_errorbar(aes(ymax = fit.1, ymin = fit.2), width = 0) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  theme_bw() +
  labs(x = "Number of groups in grid cell", y = "Coefficient Estimate") +
  theme_bw() +
  theme(panel.grid = element_blank(),
        legend.position = "none")

# * Save Plot -------------------------------------------------------------

ggsave(num_groups_plot, filename = "./plots/figure6.pdf", 
       width = 5, height = 3, scale = .9)
